### Create 9 Panel plots of FragStats with all 202 species overlaid

### Out directory

out.dir = '/home1/99/jc152199/FragStats/ImagesTBBB/'

### Directory where raw .csvs are located

csv.dir = '/home1/99/jc152199/FragStats/RawCsv2/'

### List of files in csv.dir

csv.files = list.files(csv.dir,full.names=T)

### Names of species

#spp.names = substr(list.files(csv.dir),1,nchar(list.files(csv.dir))-4)
spp.names = c('TBBB')

### Position tracker

name.track = 1

### Start looping through species

### Blank object to bind data onto

mean.2plot = NULL

for (csv in csv.files)

	{
	
	### Read in a .csv file for a single species
	
	tfile = read.csv(paste(csv,sep=''))
	
	### Save the 1990 data as a single object
	
	t1990_realized = tfile[369,c(1,6,8,10,11,16,18,23,36,42,45)]
	t1990_realized_nosmall = tfile[553,c(1,6,8,10,11,16,18,23,36,42,45)]
	
	### Remove unncessary columns
	
	tfile = tfile[,c(1,2,3,4,5,6,8,10,11,16,18,23,36,42,45)]
	
	### Omit NA's from tfile
	
	tfile = na.omit(tfile)
	
	### Subset the .csv file to include only information about the SRESA2 scenario
	
	tfile = tfile[grep('realized',tfile$dist.type),]
	
	### Aggregate realizations across GCM's to get a per year summary for each GCM
	
	t.ag = aggregate(tfile[c(7:15)],by=list(GCM=tfile$GCM,year=tfile$year,dist.type =tfile$dist.type),FUN=mean)
	
	### Create a vector the same length as nrow(t.ag) of species name
	
	t.names = rep(spp.names[name.track],nrow(t.ag))
	
	### Bind to object t.ag
	
	t.ag = data.frame(spp=t.names,t.ag)
	
	### Convert total area and mean patch area to a proportion of their 1990 size
	
	t.ag$total.area[which(t.ag$dist.type=='realized')] = t.ag$total.area[which(t.ag$dist.type=='realized')] / t1990_realized$total.area
	t.ag$mean.patch.area[which(t.ag$dist.type=='realized')] = t.ag$mean.patch.area[which(t.ag$dist.type=='realized')] / t1990_realized$mean.patch.area
	t.ag$total.area[which(t.ag$dist.type!='realized')] = t.ag$total.area[which(t.ag$dist.type!='realized')] / t1990_realized_nosmall$total.area
	t.ag$mean.patch.area[which(t.ag$dist.type!='realized')] = t.ag$mean.patch.area[which(t.ag$dist.type!='realized')] / t1990_realized_nosmall$mean.patch.area
	
	### Convert 1990 value of total area and mean patch area to 1
	
	t1990_realized$total.area = 1
	t1990_realized$mean.patch.area = 1
	t1990_realized_nosmall$total.area = 1
	t1990_realized_nosmall$mean.patch.area = 1
	
	### Append species name to 1990 data
	
	t1990_realized = data.frame(spp=spp.names[name.track],t1990_realized)
	t1990_realized_nosmall = data.frame(spp=spp.names[name.track],t1990_realized_nosmall)
		
	### Aggregate GCM's across years to calculate a yearly mean and SD for each statistic
	
	tt.ag.mean = aggregate(t.ag[5:13],by=list(spp=t.ag$spp,year=t.ag$year, dist.type=t.ag$dist.type),FUN=mean)
	
	tt.ag.sd = aggregate(t.ag[5:13],by=list(spp=t.ag$spp,year=t.ag$year, dist.type=t.ag$dist.type),FUN=sd)
	
	
	### Append the 1990 data to the mean
	
	tt.ag.mean = rbind(t1990_realized,t1990_realized_nosmall,tt.ag.mean)
	
	### Bind data for all species into one dataframe
	
	mean.2plot = rbind(mean.2plot,tt.ag.mean)
	
	### Track position
	
	cat(name.track,'\n')
	
	name.track = name.track+1
	
	years = unique(tt.ag.sd$year)
	
	### Aseemble the polygon object (a 2 column dataframe) for plotting
	
	poly.out = NULL
	
	for (dist.type in c('realized','realized.NO.small.patches'))
	
		{
	
	for (yearx in years)
		
		{
		
		ttt.mean = tt.ag.mean[which(tt.ag.mean$year==yearx & tt.ag.mean$dist.type==dist.type),]
		
		ttt.sd = tt.ag.sd[which(tt.ag.sd$year==yearx & tt.ag.sd$dist.type==dist.type),]
		
		t.year.min = data.frame(dist.type=dist.type,year=yearx,ttt.mean[,c(4:12)] - 1.96*ttt.sd[,c(4:12)])
		
		if (yearx==2000 & dist.type=='realized') {poly.out = rbind(poly.out,t1990_realized[2:12],t.year.min)}
		
		if (yearx==2000 & dist.type=='realized.NO.small.patches') {poly.out = rbind(poly.out,t1990_realized_nosmall[2:12],t.year.min)}
		
		if (yearx!=2000) {poly.out = rbind(poly.out,t.year.min)}
		
		}
		
	for (yearx in rev(years))
		
		{
		
		ttt.mean = tt.ag.mean[which(tt.ag.mean$year==yearx & tt.ag.mean$dist.type==dist.type),]
		
		ttt.sd = tt.ag.sd[which(tt.ag.sd$year==yearx & tt.ag.sd$dist.type==dist.type),]
		
		t.year.max = data.frame(dist.type=dist.type,year=yearx,ttt.mean[,c(4:12)] + 1.96*ttt.sd[,c(4:12)])
		
		poly.out = rbind(poly.out,t.year.max)
		
		
		
		}
	
		}
	
### Close loop and begin plotting
	
### Proper names for plots
	
stat.names = c('Abundance','Patches','Total Area','Landscape Shape Index','Mean Patch Area','Mean Perimeter Area Ratio','Landscape Core','Aggregation Index','Eff. Mesh Size')

### Edit names for mean.2plot

names(mean.2plot)[4:12]=stat.names
names(poly.out)[3:11]=stat.names
	
### Proper names for y-labels
	
ylab.names = c('Proportion','Count','Proportion','Percent','Proportion','Mean Ratio * 10^4','Proportion','Percent','Continuous / 10^6')

#######################################
### Define species codes for species of interest as a vector

sois = c('TBBB')

### Calculate y-lims for all stats using a loop

### Deine ylims as a NULL object

ylims = NULL

for (colx in 3:11)

	{
	
	t.data = data.frame(poly.out[,1:2],poly.out[,colx])
	
	#t.data = t.data[which(t.data$spp %in% sois),]
	
	t.min = min(t.data[,3])
	
	t.max = max(t.data[,3])
	
	t.out = data.frame(stat=stat.names[colx-2],min=t.min,max=t.max)
	
	ylims = rbind(ylims,t.out)
	
	}
	
### Ylims dataframe assembled

### Plot

### First subset by frag.stat

### Position tracker
		
plot.track=1
	
	for (stat in stat.names)
		
		{
		
		### Start the .png device driver
		
		png(paste(out.dir,stat,' TBBB.png',sep=''), units='in', width=10, height=10, res=300)

		### Ylims for a single stat
		
		t.ylims = ylims[which(ylims$stat==stat),]
		
		### Mean data for a single stat
		
		t.2plot = data.frame(dist.type=mean.2plot$dist.type,year=mean.2plot$year, stat=mean.2plot[,which(names(mean.2plot)==stat)])
		
		# Poly data for a single stat
		
		t.poly = data.frame(dist.type = poly.out$dist.type,year=poly.out$year, stat=poly.out[,which(names(poly.out)==stat)])
		
		### Create the axis tick positions and labels for each plot (conditional on stat type)
		
			if (stat=='Abundance') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,2)}
			if (stat=='Patches') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,0)}
			if (stat=='Total Area') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,2)}
			if (stat=='Landscape Shape Index') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,0)}
			if (stat=='Mean Patch Area') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,2)}
			if (stat=='Mean Perimeter Area Ratio') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4)*1000;t.yticks = round(t.yticks,0);t.2plot$stat=t.2plot$stat*1000;t.poly$stat=t.poly$stat*1000;t.ylims[1,2:3]=t.ylims[1,2:3]*1000}
			if (stat=='Landscape Core') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,2)}
			if (stat=='Aggregation Index') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4);t.yticks = round(t.yticks,2)}
			if (stat=='Eff. Mesh Size') {t.yticks = seq(t.ylims[1,2],t.ylims[1,3],(t.ylims[1,3]-t.ylims[1,2])/4)/1000000;t.yticks = round(t.yticks,0);t.2plot$stat=t.2plot$stat/1000000;t.poly$stat=t.poly$stat/1000000;t.ylims[1,2:3]=t.ylims[1,2:3]/1000000}

	
		### Configure the plot space
		
		plot(t.2plot[,2],t.2plot[,3],xlim=c(1990,2080),ylim=c(t.ylims[1,2],t.ylims[1,3]),main=paste(stat),type='n',axes=F,xlab='Year',ylab=ylab.names[plot.track])
		
		### Color object
		
		col.ob  =c('#0000FF','#00FF00')
		
		### Position tracker
		
		col.track = 1
		
			for (dist.type in c('realized','realized.NO.small.patches'))
	
				{
	
				### Subset to data for one one patch size only
				
				tt.2plot = t.2plot[which(t.2plot$dist.type==dist.type),]
				
				### Subset to data for one patch size only
				
				tt.poly = t.poly[which(t.poly$dist.type==dist.type),]
				
				### Overlay the mean of 1 statistic for 1 species
		
				lines(tt.2plot$year,tt.2plot$stat, col=col.ob[col.track],lwd=1)
				
				### Poly the polygon shape for 1 stat
				
				polygon(tt.poly[,2],tt.poly[,3], col=paste(col.ob[col.track],25,sep=''), lty=1, lwd=1, border='lightblue')
				
				### Legend for plot
		
				legend('topright',bty='n',legend=c('Realized All Patches','Realized Excluding Small Patches'),text.col=c('#0000FF','#00FF00'))
				
				### X-axis
				
				axis(side=1,at=c(1990,2000,2010,2020,2030,2040,2050,2060,2070,2080),labels=c(1990,'',2010,'',2030,'',2050,'',2070,''))
				
				### Y-axis
				
				axis(side=2,at=c(t.yticks),labels=c(t.yticks))
				
				### Position trackers
				
				col.track = col.track+1
				
				### Report progress
				
				#cat(spp,' - ',stat,'\n',sep='')
				
				}
		
			### All 9 stats plotted for a single species
		
			cat(stat,'\n')
			
			dev.off()
			
			plot.track=plot.track+1
		
	}
	
}
	
# Done
		
		
	
	